Mapping in R

Mapping

This workshop material is a slightly modified varient of the mapping section from the SCAR-EGABI Tools for Southern Ocean Spatial Analysis and Modelling workshop written by Ben Raymond. The full workshop material is available here: https://scar.github.io/EGABIcourse19/index.html

Here we are going to explore a range of ways to create maps; from the original mapping packages in R (maps), making simple maps in ggplot2, and then creating polar maps using SOmap.

The first thing to know about mapping in R is perfectly summarised by @djnavarro:

Maps in R

The oldest and most core general mapping package in R is the maps package. It has a simple whole-world coastline data set for immediate use.

The data underlying this live map is available by capturing the output as an actual object. Notice that the coastline for Antarctica does not extend to the south pole, and that parts of Russia that are east of 180 longitude are not in the western part of the map.

A very similar and more modern data set is available in the maptools package.

This data set aligns exactly to the conventional -180/180 -90/90 extent of the longitude/latitude projection.

Exercises

  1. How can we find the longitude and latitude ranges of the maps data m and the maptools data wrld_simpl?
    Show answer

    Answer 1: range(m$x, na.rm = TRUE) range(m$y, na.rm = TRUE) also m$range

  2. Can we draw polygons with a fill colour with the maps package?
    Show answer

    Answer 2: polygon(lonlat, col = "grey") does not work, and map(mp, fill = TRUE, col = "grey") does not work, but maps::map(fill = TRUE, col = "grey") does seem to work.

What’s going on? Look at the very south-eastern corner of the map. The “coastline” has been extended to the very south boundary of the available area.

When we add the old maps coastline see that it does not extend to 90S and it does not traverse the southern boundary.

One reason for this is that if we choose a projection where the east and west edges of the Antarctic coastline meet then we get what looks a fairly clean join.

If we try the same with wrld_simpl it’s not as neat. We have a strange “seam” that points exactly to the south pole (our projection is centred on longitude = 0, and latitude = -90.

Let’s use the maps data!

In m we have the maps data structure, and this looks promising.

## List of 4
##  $ x    : num [1:82403] -69.9 -69.9 -69.9 -70 -70.1 ...
##  $ y    : num [1:82403] 12.5 12.4 12.4 12.5 12.5 ...
##  $ range: num [1:4] -180 190.3 -85.2 83.6
##  $ names: chr [1:1627] "Aruba" "Afghanistan" "Angola" "Angola:Cabinda" ...
##  - attr(*, "class")= chr "map"
## [1] -12709814  12704237 -12576156  12470787

The problem is that the maps database has enough internal structure to join lines correctly, with NA gaps between different connected linestrings, but not enough to draw these things as polygons. A similar problem occurs in the default projection. While wrld_simpl has been extend by placing two dummy coordinates at the east and west versions of the south pole, this data set does not have those.

We have to look quite carefully to understand what is happening, but this is wrapping around overlapping itself and so close to the southern bound we barely notice.

## [1] -20037508  20037508 -20179524  18351859

ggplot2

So far all the maps we have made use graphics::plot() but we can also use ggplot2 to make our maps.

As a quick example of how we might do this:

##         lon      lat group    id
## 1 -69.89912 12.45200     1 Aruba
## 2 -69.89571 12.42300     1 Aruba
## 3 -69.94219 12.43853     1 Aruba
## 4 -70.00415 12.50049     1 Aruba
## 5 -70.06612 12.54697     1 Aruba
## 6 -70.05088 12.59707     1 Aruba

The map_data function is an internal ggplot2 function to access the data in the maps package. Above we create a world object with all the data from the countries of the world.

Below we plot this using ggplot2 plotting the countries as polygons. Note the use of group to specify the breaks in the data to make different countries different polygons.

Lonlat

Stereographic

Orthogonal

Adding some data? We probably need some first so we will make some fake data quickly

To add the points we use geom_point() on the end of our map.

SOmap

The SOmap package is intended to solve some of these problems, and provide an easier way to produce nice-looking maps of Antarctica and the Southern Ocean. It is primarily focused on maps in polar stereographic projection (although the SOmap_auto function extends this to other projections). SOmap won’t necessarily get you exactly the map you want in every circumstance, but the idea is that in most cases it should get you close enough, and if need be you can make modifications to suit your exact purposes.

Please bear in mind that SOmap is still in development, and so its functionality (function parameters and/or behaviour) may change.

By default, SOmap works with base graphics (and associated functionality from packages such as raster and sp). It is also possible to work with ggplot2-based graphics, as described below.

Start by installing the SOmap package if you haven’t done so already:

remotes::install_github("AustralianAntarcticDivision/SOmap")

Temporary note: the code here was written using the development version of SOmap. You may need to install the development version using remotes::install_github("AustralianAntarcticDivision/SOmap", ref = "dev-0.5").

Then load the package:

## Loading required package: raster
## Loading required package: sp

Circumpolar maps

A basic circumpolar map in polar stereographic projection:

## Loading required namespace: rgeos

SOmanagement() provides a number of contextual layers such as MPA boundaries and management zones.

Adding raster layers

First let’s construct some artificial raster data (in longitude-latitude space) for demonstration purposes:

SOplot will reproject and plot this for us:

The legend is out of character with the rest of the map. We can use SOleg to fix that:

OK, well that worked but clearly the labels need tidying up. The easiest way is probably to set the number of decimal places in the label values via the rnd parameter:

Alternatively, we could explicitly set the colour range and labels.

Note that if we don’t want to show the bathymetric legend, we may run into problems:

The legend has been chopped off because the layout has not left enough space around the map for the curved legend. Currently, the best solution is probably to generate the SOmap object with the bathymetric legend, but then remove it before plotting (see the Modifying map objects section for more details on this):

Multiple rasters:

####Create a raster from points

A common situation is recieving a large amount of overlapping point data and wanting to look at it. One option is to plot the points on the map:

## No projection provided, assuming longlat

This however is not particularly useful if you want to know much more than the bounds of points. If we want to know how often an area is visited than we want to use a raster. There is a (currently) internal function in SOmap called SObin which will take your lat long points, reproject them and create a raster layer from them.

## No projection provided, assuming longlat

Non-circumpolar maps

The SOmap_auto function will take your input data and make a guess at an appropriate projection and extent to use. Note that this is not always going to guess the best projection and extent, so you should view it as a starting point from which you can generate a map to your exact requirements.

Use the elephant seal track data bundled with the package:

Just a blank map to which you could add other things:

You can pass a raster as input data, but note that it won’t plot the raster (it uses its extent to infer an appropriate extent for the map):

But we can add the raster if we wish:

We can force a particular projection:

Same but by supplying a full proj4 string to target:

See the SOmap_auto vignette for more examples.

Plotting via ggplot2

The SOmap and SOmap_auto functions do their plotting using base graphics. If you are more comfortable working with ggplot2, this is also possible. The SOgg function takes an object created by one of those functions (using base graphics) and converts it to use ggplot2 graphics instead. As with other SOmap functions, this returns an object (of class SOmap_gg or SOmap_auto_gg) that contains all of the information needed to generate the map. Printing or plotting this object will cause it to construct a ggplot object. Printing or plotting that object will cause it to be drawn to the graphics device, just like any other ggplot object.

## [1] "SOmap_gg"
## [1] "gg"     "ggplot"

Or in one step (this will cause myplot to be converted to SOmap’s internal gg format, then a ggplot object constructed from that, then that object will be plotted):

Modifying map objects (advanced usage)

The goal of SOmap is to make it fairly easy to produce a fairly-good-looking map that will be adequate for most mapping requirements. It will never be possible to automatically produce a perfect map in every circumstance, but the aim is to have a low-effort way of getting fairly close most of the time.

This section describes some approaches to modifying a map to get it closer to your particular needs. Be warned: depending on the exact modifications needed, this might get you pretty close to the crumbling edge of SOmap development. In particular, anything that requires modifying the internal structure of an SOmap object may change in the future (with luck, we’ll make this sort of thing easier - but we’re not there yet.)

Modifying base graphics maps

Calls to SOmap(), SOmanagement(), SOmap_auto() return an object of class SOmap, SOmap_management, or SOmap_auto. These objects contain all of the data and plotting instructions required to draw the map. Calling print() or plot() on one of these objects will cause that code to be executed, and the object to be drawn in the current graphics device. Hence, calling SOmap() directly without assigning the result to a variable will make it appear in the graphics device, because the returned object is being printed to the console (and thus triggering the print method). But you can also assign the result to a variable, e.g. myplot <- SOmap() and then explicitly plot the object with plot(myplot). The advantage of this is that you can potentially manipulate the myplot object to make changes to the map before plotting it.

Note, this is likely to be fragile. Proceed at your own risk!

##  [1] "projection"    "target"        "straight"      "trim"         
##  [5] "bathy"         "box"           "plot_sequence" "coastline"    
##  [9] "ice"           "outer_mask"    "bathy_legend"  "border"

The object contains a plot_sequence component, which defines the order in which each part of the plot is drawn. The other components of the object contain the code required to draw each part. Take e.g. the ice component (this is the ice shelves, glacier tongues, etc). This is a list (in this case with only one element). Each element of the list specifies a function to run along with arguments to pass to it:

## List of 1
##  $ :List of 2
##   ..$ plotfun : chr "plot"
##   ..$ plotargs:List of 4
##   .. ..$ x     :sfc_POLYGON of length 354; first list element: List of 1
##   .. .. ..$ : num [1:5, 1:2] 1022981 1026000 1021994 1021935 1022981 ...
##   .. .. ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
##   .. ..$ col   : logi NA
##   .. ..$ border: chr "black"
##   .. ..$ add   : logi TRUE
##   ..- attr(*, "class")= chr "SO_plotter"

We can modify the function and/or its arguments:

We can remove entire components:

But note that some elements are required. In particular, the bathymetry layer can’t currently be removed because the code that draws this is also the code that creates the plot page via plot.new(). The code below would fail outright if there was no existing existing plot. If there was an existing plot in the graphics device, this code would run but give unpredictable results because it would draw on top of the previously-setup plot:

One way around this would be to simply replace all of the bathymetric data values with NAs. The plotting code would still have the extent of the bathymetric layer that it needs in order to set up the plot, but no data would be shown:

## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf

We could also replace the bathymetry data with another raster object. Note that we do need to be careful about the extent and projection of this raster. For example, replacing the bathymetry raster with the ice raster (which has the same polar stereographic projection but smaller extent) gives:

It’s chopped off because the extent of the ice raster is being used to set the plot extent. But if we extend the ice raster to match the map extent:

Modifying ggplot maps

We can modify ggplot2-based maps at two levels.

Modifying the ggplot object.

Remember that printing or plotting a SOmap_gg object produces a ggplot object. This can be modified by adding e.g. layers or themes just like a normal ggplot. Remember to load the ggplot2 library now that we are using ggplot2 functions directly.

Multiple rasters or multiple sets of points gets tricky if they are on different scales, because ggplot2 is only designed to work with a single colour scale per geometry type. You can try your luck with the ggnewscale or relayer packages, although both are in a fairly experimental stage of development.

Modifying the SOmap_gg object

SOmap_gg objects are similar in structure to SOmap objects, in that they contain all of the data and plotting instructions required to draw the map:

##  [1] "projection"    "target"        "straight"      "trim"         
##  [5] "init"          "bathy"         "coord"         "plot_sequence"
##  [9] "scale_fill"    "bathy_legend"  "coastline"     "ice"          
## [13] "axis_labels"   "theme"         "border"

However, instead of base plotting functions, SOmap_gg objects use ggplot2 function calls, e.g.:

## [1] "ggplot2::geom_sf"

We can modify these functions and/or arguments in a similar manner to SOmap objects.

Or remove the bathymetric raster layer:

Or replace it with a different raster (use the ice raster as an example):

Supporting data for maps

When constructing maps, we commonly want to show features like oceanographic fronts, ice extent, coastline, place names, and MPA boundaries. There are a few sources of such data:

  • some layers are bundled into SOmap, see the SOmap::SOmap_data object
  • antanym provides access to the SCAR Composite Gazetteer of place names
  • the quantarcticR package provides access to Quantarctica data layers.

quantarcticR

Note, this package is still in development, so the usage as shown here might change in later versions. Install if needed:

Example usage:

## # A tibble: 6 x 5
##   layername          main_file                       type   cached download_size
##   <chr>              <chr>                           <chr>  <lgl>    <fs::bytes>
## 1 Overview place na~ "C:\\Users\\dale_mas\\AppData\~ shape~ FALSE         19.74K
## 2 COMNAP listed fac~ "C:\\Users\\dale_mas\\AppData\~ shape~ FALSE        691.92K
## 3 Subantarctic stat~ "C:\\Users\\dale_mas\\AppData\~ shape~ FALSE        691.92K
## 4 SCAR Composite ga~ "C:\\Users\\dale_mas\\AppData\~ shape~ FALSE        329.05M
## 5 IBO-IOC GEBCO Fea~ "C:\\Users\\dale_mas\\AppData\~ shape~ FALSE          1.25M
## 6 IBO-IOC GEBCO Fea~ "C:\\Users\\dale_mas\\AppData\~ shape~ FALSE          1.25M
## # A tibble: 1 x 11
##   layername datasource layer_attributes srs_attributes provider abstract extent
##   <chr>     <chr>      <list>           <list>         <chr>    <chr>    <list>
## 1 Median s~ SeaIce/Me~ <NULL>           <tibble [1 x ~ ogr      "Monthl~ <NULL>
## # ... with 4 more variables: palette <list>, type <chr>,
## #   download_size <fs::bytes>, main_file <chr>

Citations

No matter what packages you use in R you should always cite them. Often packages will have far more time invested in them than some publications and the authors deserve acknowledgement. R makes this easy for example, the citation for SOmap can be found with:

## Loading required namespace: rgeos
Fig. inf: shiny Southern Ocean map made with SOmap (Maschette et al., 2019).

Fig. inf: shiny Southern Ocean map made with SOmap (Maschette et al., 2019). 

## 
##   Maschette, D., Sumner, M., Raymond, B. (2019). SOmap: Southern Ocean
##   maps. Version 0.5.0
##   https://github.com/AustralianAntarcticDivision/SOmap
## 
## A BibTeX entry for LaTeX users is
## 
##   @Manual{,
##     title = {SOmap: Southern Ocean maps},
##     author = {Dale Maschette and Michael Sumner and Ben Raymond},
##     year = {2019},
##     url = {https://github.com/AustralianAntarcticDivision/SOmap},
##     version = {R package version 0.5.0.9999},
##   }
## 
## To cite ggplot2 in publications, please use:
## 
##   H. Wickham. ggplot2: Elegant Graphics for Data Analysis.
##   Springer-Verlag New York, 2016.
## 
## A BibTeX entry for LaTeX users is
## 
##   @Book{,
##     author = {Hadley Wickham},
##     title = {ggplot2: Elegant Graphics for Data Analysis},
##     publisher = {Springer-Verlag New York},
##     year = {2016},
##     isbn = {978-3-319-24277-4},
##     url = {https://ggplot2.tidyverse.org},
##   }

Dale Maschette, Ben Raymond & Michael Sumner

02/10/2019